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Abstract. The two dimensionality plus the linear band structure of graphene 
leads to new behavior of the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction, 
which is the interaction between two magnetic moments mediated by the electrons 
of the host crystal. We study this interaction from linear response theory. There 
are two equivalent methods both of which may be used for the calculation of the 
susceptibility, one involving the integral over a product of two Green's functions 
and the second that involves the excitations between occupied and unoccupied 
states, which was followed in the original work of Ruderman and Kittel. Unlike 
the J oc (2kpR)~ 2 sin(2kpR) behavior of an ordinary two-dimensional (2D) metal, 
J in graphene falls off as 1/i? 3 , shows the 1 + cos((K — K').R)-type of behavior, 
which contains an interference term between the two Dirac cones, and it oscillates 
for certain directions and not for others. Quite interestingly, irrespective of any 
oscillations, the RKKY interaction in graphene is always ferromagnetic for moments 
located on the same sublattice and antiferromagnetic for moments on the opposite 
sublattices, a result that follows from particle-hole symmetry. 
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I. INTRODUCTION 

Graphene, which is a plane of carbon atoms with a honeycomb lattice, is of considerable 
interest 1 owing to its two- dimensionality and a linear band structure as opposed to the 
quadratic band structure in typical materials. These features also introduce new behaviors in 
the RKKY interaction, which has been extensively studied beginning with the original works 
of Ruderman and Kittel 2 , Kasuya 3 , and Yosida. 4 The RKKY interaction is the interaction 




FIG. 1: (Color online) The RKKY interaction between two magnetic moments mediated by the host 
electrons in the crystal. 

between two magnetic moments mediated by the conduction electrons in the host material. 
The first moment perturbs the conduction electrons, which is seen by the second moment 
leading to an indirect exchange interaction as illustrated in Fig. (1). For the free electron 
gas with quadratic bands E = h 2 k 2 /2m and Fermi momentum hp, the strength of the 
interaction, J, is given by the expression 2-6 



where Si(x) is the sine integral function and x = 2k pR with R being the separation distance 
between the two moments. As can be seen from Eq. (1), the common behavior of J in any 
dimension is characterized by the power-law decay with some oscillations whose period is 
scaled by the Fermi momentum kp. 

In contrast to these simple forms, the form of the RKKY interaction in graphene is quite 
complex and has been the subject of many papers both for the doped and undoped cases. 7-15 
Owing to the lattice structure and a gapless density of states at the Fermi energy with the 
linear bands occurring at two different points (the Dirac points) in the Brillouin zone (BZ), 
the RKKY interaction depends on the directionality as well as on the sublattice locations of 
the two magnetic moments and, in addition, it contains an interference term coming from 
the two Dirac points K and K'. As we show in this article, the net result is 
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where C is a constant, 6r is a position-dependent phase angle, and the subscripts in J a p 
indicate the sublattice location of the moments. The detailed method of how to obtain these 
results was presented in our earlier study, where we used an expression for the susceptibility 
in terms of an integral over the product of two Green's functions. In that method, a cut-off 
function 7,11 was necessary to evaluate the integrals, although quite satisfactorily the final 
results did not depend on the exact cut-off function. 

It is illuminating to evaluate the RKKY interaction using an alternative expression for 
the response, expressed in terms of the excitations of the system, a method familiar in the 
literature from the original Ruderman-Kittel formulation 2 and one that has been recently 
applied to graphene as well. 15 It yields the same results we had obtained before, without 
necessitating the use of a cut-off function. In addition to the details of the method, we 
discuss the salient features of the results for the RKKY interaction in graphene. 



The RKKY interaction is directly proportional to the susceptibility. The response of the 
charge density, n, to a perturbing potential, V, may be written in terms of the integral over 
the unperturbed Green's function 



where Xap{ r i r ') = Sn a (r) / '8Vp(r') is the charge susceptibility for a crystal with the Greek 
subscripts indicating the sublattice indices, 6Vp(r') is a spin-independent perturbing poten- 
tial and 5n a (r) is the induced charge density that includes both spin channels. 

We outline briefly the derivation of the alternative expression for the susceptibility in 
terms of the energy excitations. This may be obtained by using the spectral representation 
of the Green's function 



where ip^ s is the sublattice component of the unperturbed eigenfunction with the corre- 
sponding energy For a crystalline structure, {k, s} denotes the Bloch momentum and 
the band index; else, it just denotes a complete set of states. Plugging Eq. (5) into the 
expression Eq. (4), one finds the result 



II. SUSCEPTIBILITY 
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It can be easily shown that under the interchange of ks and k's', the real part of the product 
of the four wave functions appearing in the equation above is even, while its imaginary part 
is odd, and at the same time, both the real and the imaginary parts of the product of the 
momentum-space Green's function are even. This makes the second line zero. In addition, 
in order to produce a final compact equation, we replace the real part in the first line by the 
entire complex quantity, as the extra term introduced thereby gives a zero net result when 
summed. We therefore obtain the expression 

Xa p(r, r') = Y,^s{rWl{r')^Ar^{r) x(ks, k's'), (7) 
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The last equality is obtained by using the relationship lim ?? _ ) . +(:r ±i?7) -1 = V(l/x) ^fiirS^x). 
Clearly, the integral is non-zero only if the eigenstate ks is occupied while k's' is empty 
or vice versa. Thus, corresponding to these two processes, Eqs. (7) and (8) lead to two 
terms, which may be combined into a compact expression by using the Fermi function 
f(e) = 9(e F — e), where the step function 9(x) is, as usual, 1 if x > and otherwise. This 
leads to our desired result 

Xa ,(r,r') = 2 ^/!^L^) f b (r)tf (04(0^(r)- (9) 
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This is a well-known formula in the linear response theory and is the central equation in 
this paper, from which we will compute the sublattice susceptibilities for graphene. Note 
that under the interchange of ks and k's', the summand in Eq. (9) goes into its complex 
conjugate, so that only the real part survives in the summation. 



As already stated, the RKKY interaction can be expressed in terms of the susceptibility. 
Taking the interaction between the localized moments and the conduction electrons as a 
contact interaction in the form 

V=-\ (S 1 -s 1 + S 2 -s 2 ), (10) 

where Sj is the conduction electron spin density at site i, the interaction energy between the 
two localized moments may be written as 16 

E(r,r') = J af} (r,r , )S 1 -S 2 , (11) 

where r, r' denote the lattice positions of the two spins and the RKKY interaction J a p(r, r') 
is simply proportional to the susceptibility 
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J a p{r,r') = —Xap{ry). (12) 



III. THE RKKY PROBLEM IN GRAPHENE 



We consider the nearest-neighbor tight-binding Hamiltonian for the ir electrons in 
graphene with the interaction between the host electrons and the localized magnetic mo- 
ments given by Eq. (10) as before. Thus we have 



= -t^24 a c ja + H.c, 
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in which (ij) denotes summation over distinct pairs of nearest-neighbor sites with hopping 
parameter t and cj CT (c^) is the creation (annihilation) operator for an electron with spin 
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FIG. 2: (Color online) The graphene honeycomb lattice with two different sublattices, shown as red and 
blue dots. The figure also shows the corresponding BZ with the Dirac points K and K' and two common 
directions in the direct lattice (zigzag and armchair). Ti and T 2 are the two primitive translation vectors 
of the direct lattice and the three nearest-neighbor distance vectors are indicated by d\, d 2 , and d 3 . 



index a and combined site-sublattice index i. Two set of Bloch sums in the momentum- 
sublattice representation, viz., \ka) = c^jO) are introduced to diagonalize T-Lq. The Block 
sums are 
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where N is the number of unit cells in the lattice, Ri denotes the cell positions, and r a 
denotes the basis atom positions in the unit cell. We take ta = and tb = d 1 . In the basis 
of the sublattice Bloch wave functions \ka), the Hamiltonian becomes 

f(k) 
J*(k) 

where f{k) = —t X]j=i e * fc dj ' with dj being the nearest-neighbor position vectors and a 
being the carbon-carbon bond length. The unperturbed eigenstates of Hk are given by 

£ks = s\f(k)\ 
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where the band index s = ± denotes the conduction and the valence bands and the phase 
factor <f)(k) = e l ® k = f(k)/\f(k)\. Note that by choosing the Bloch sum (14) that includes 
the phase factor e ik ' Ta , the basis wave function for the two atoms in a particular cell also 
contains the phase factors e tk ' Ta . 

Many of the fascinating properties of graphene like high electron mobility are attributed 
to its cone-shaped, linear band structure around the corners of the BZ, called the Dirac 




FIG. 3: (Color online) The direction of pseudo-spin around Dirac points K and K' . The wave functions have 
definite helicities ±1 depending on whether the pseudo-spin er is parallel or antiparallel to the momentum 

q 



cones. The expansion of function f(k = q + K D ) for small q around all six Dirac points 
takes the form f(q + Kp) = v F q <fi(q), leading to the linear band structure 

e = ±v F q (17) 

with Fermi velocity v F = 3ta/2. The phase factors appearing in the wave function Eq. 
(16) are, however, different near different Dirac points. For all six Dirac points shown 
in Fig. (2), starting from the top points and going counter-clockwise, these phase are: 

0(g) = e^/ 3 -^), -e^+f^ _ e -*? 9) e i(2n/3+8 g )^ _ e i(2n/3-8 g ) ? e iO q with pol&r &ngle Q f g 

defined as 6 q = tan -1 (q y / q x ) ? Note that these phases are a direct consequence of using the 
extra phase factor e tk - Ta in the Bloch sum, Eq. (14). This choice is preferred as all physical 
quantities will be evaluated at the actual positions of the atoms rather than the unit cell 
positions Ri in which a particular atom is located. However, the second choice of Bloch 
sums, not adopted in this paper, will have the same phase 0(g) near all Dirac points K or 
K' and is used sometimes in the literature. 

Near the Dirac cones, the Hamiltonian assumes a simple form 

K q +K = v F (T* ■ q, U q+K > = -v F a ■ q, (18) 

leading to its interpretation in terms of the pseudo-spins. Here cr = (a x , a y ) are the pseudo- 
spin Pauli matrices describing the two sublattices and cr* = (a x , —cr y ). The two-component 
central-cell wave functions are 
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in the basis set of the carbon orbitals without any phase factors e lk ' Ta in their definitions. The 
wave functions have definite helicities ±1 corresponding to the eigenvalues of the operator 
h = <r ■ q as indicated in Fig (3). 



A. Moments on the Same Sublattice: Jaa{R) 

We will now use these eigenstates to evaluate the RKKY interaction using the suscepti- 
bility expression Eq. (9). Using Eq. (16) in Eq. (9) for both moments on A-sublattice, one 
at the origin and the other at the atom position R, we find 

Y e -i(k-k')-R 

k,k' fc lt + 

£fc-<£F<£fc' + 



We now evaluate Eq. (20) for the linear Dirac bands. We can construct the BZ such that 
it encloses two top Dirac points shown in Fig (2) and perform the summations over k and 
k' over two circles centered at the Dirac points. We first perform the k' summation, which 
yields 

1 . , pik'R i piqi-R 1 „ P iq2-R 

_Ly^_^ = e iK.R(J_y* e x +e iK'-R/±_ e , _ 

N^-?e k _-e k , + K N^e k . -v F q x > K N ^ e k _-v F q 2 >' 

k' <ji <J2 

where qi, q 2 denote the momentum with respect to the two Dirac points. Using the Jacobi- 
Anger expansion 17 for the exponentials, viz., 

oo 

e ±iq R = MqR) + 2 ^(±{) n J n (qR) cos[n(0, - 9 R )}, (22) 

71=1 

where J n (x) is the Bessel function and the integral e ±iq R d9 q = 2rrJ (qR), Eq. (21) 
yields 

Here we have used iV _1 — >■ (2f2#z) 1 / cPg with Q BZ being the area of the BZ. We now 

Q J 

perform the fc-summation in Eq. (20) using Eq. (23) following the similar steps as above 
and finally get the expression 

XAA{0 , R) = ^A ^osm lAA] (24) 

v F il BZ R 6 
Jo Jo x -\- x 

where x = qR, x' = q'R, and no cutoff has been used for the Dirac cones. To evaluate Iaa 
we define the function 

H(s) = re- s ^ XX ' Mx)Mx,) dxdx>. (26) 
Jo Jo x -\- X 

The derivative of H(s) gives the square of a Laplace transform, which can be easily evaluated, 

with the result 

^ --[fe-xMx)} 2 = -{C\xUx)]f = -j^. (27) 

Integrating this and determining the constant of integration from the condition H(0) = Iaa 
(see Eqs. (26) and (25)), we find 

H{s) = WTW-l tan ^ s + lAA - (28) 
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FIG. 4: (Color online) RKKY interaction Jaa between two moments on the same sublattice located either 
along the zigzag or the armchair directions as obtained from Eq. (29). Note that Jaa is negative for all 
R indicating a ferromagnetic interaction. Consistent with Eq. (29), J is oscillatory for the zigzag direction 
and smoothly decays along the armchair direction, while always remaining ferromagnetic. 

From the definition of H(s) in Eq. (26) we see that lim^oo H(s) = and evaluating the 
right hand side for s = oo, we immediately find Iaa = vr/16. Plugging this result for Iaa 
into Eqs. (24) and (12), we find the RKKY interaction to be 



where C = 9A 2 ft 2 /(2567rt). As C > 0, Eq. (29) represents a ferromagnetic coupling between 
the moments on the same sublattice. One can simply show that Dirac-cone oscillatory fac- 
tor, 1 + cos[(K — K') ■ R] takes the sequence of triplets of 2, 1/2, 1/2, ... with distance R 
along the zigzag direction, and becomes always 2 for the armchair direction. These results 
are shown in Fig. (4). This is consistent with the conclusion that 11 in the presence of the 
particle-hole symmetry (which is true for a bipartite lattice with no interaction between the 
members of the same sublattices), the RKKY interaction between two moments placed on 
the same sublattice is ferromagnetic, while those placed on the opposite sublattices is an- 
tiferromagnetic. Presence of the second nearest-neighbor interaction breaks this symmetry, 
which is relatively weak in graphene. 18 




(29) 



B. Moments on the Opposite Sublattices: Jab{R) 



For two moments located on the opposite sublattices, the first on the A sublattice atom 
at the origin and the second at the atom position R on the 5-sublattice, Eqs. (9) and (16) 



AB-Zigzag 



AB-Armchair 



10" 



CM 

^ 10" 5 



£ 10" 8 
10 



10 



-14 




50 100 150 200 250 
R/a 




CD 

-? 10" 7 



100 150 200 

R/a 



250 



FIG. 5: (Color online) Jab for the zigzag direction calculated from Eq. (33). Note that Jab is positive for 
all R indicating an antiferromagnetic coupling for all distances, even though the magnitude can oscillate. 
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Following similar algebra as in the previous section and performing the angle integration 
after the Jacobi-Anger expansion Eq. (22), viz., J^ 77 e ±tq ' R e ±teq d9 q = ±2iriJi(qR)e ±ldR , the 
susceptibility is given by 
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where x = qR and x' = q'R. This integral can be evaluated following the above method 
of Laplace transform with replacing Jq(x) by J\{x) in Eqs. (26) and (27). The result is 
Iab = 37r/16. With this, Eqs. (31) and (12) yield the RKKY interaction 

1 + cos[(K -K')-R + tt- 26 R ] 



Jab(R) = 3C x 
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Clearly, Jab is always antiferromagnetic as required by particle-hole symmetry, even though 
its magnitude may oscillate with distance. The results are plotted in Fig. 5 for two different 
directions in the graphene lattice. 



IV. SUMMARY 



We have studied the RKKY interaction in graphene by computing the susceptibility Eq. 
(9) in terms of the one-particle excitations in the system following the original method of 



Ruderman and Kittel. 2 The results are the same as evaluated using the integration over the 
product of the Green's function Eq. (4), a method we have adopted in our earlier works. 7,8 
The present method is somewhat better for graphene in that no cut-off functions are needed 
to perform the integrals. 

The RKKY interaction in graphene has several interesting features. Unlike the 1/R 2 
fall-off of the RKKY interaction in the standard 2D metals with quadratic dispersion, the 
linear band structure of graphene leads to the 1/R 3 dependence on distance. The inter- 
ference between the two Dirac cones in the Brillouin zone produces an interference effect 
that leads to an oscillatory behavior of the RKKY interaction. However, even though the 
magnitude may oscillate with distance along certain directions, the sign of the interaction 
is always ferromagnetic for moments located on the same sublattice and antiferromagnetic 
for moments located on the opposite sublattices, a result that follows from the particle-hole 
symmetry of graphene. 11 Lately it has been possible to electron or hole dope graphene by a 
gate voltage. The same analysis can be extended to this case. Some results for the RKKY 
interaction for the doped system have been presented in the literature. 8 

We finally note that there is currently considerable interest in the topological insulators, 
where one also has two-dimensionality and a linear band structure. However, there are im- 
portant differences between the topological insulators and graphene. For instance, graphene 
contains an even number of Dirac cones (four including spin and valley degeneracies), while 
the topological insulators contain an odd number of cones. Secondly, in graphene, the linear 
band structure originates by the presence of a pseudo-spin, while in the topological insula- 
tors we have real spins, so that a magnetic impurity opens up a local gap and suppresses 
the local density of states. As in the usual Fermi liquid, if the surface state has a finite 
Fermi wave vector fop, the sign of the RKKY interaction oscillates with wavelength oc /jp/2. 
However, if the Fermi level is close to the Dirac point, the RKKY interaction will always be 
ferromagnetic as a uniform spin polarization can maximize the gap opened on the surface. 19 
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